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f^ ' We study statistical properties ol spatial distances between successive earthquakes, the so-called 

^^ ' hypocenter intervals, produced by a two-dimensional (2D) Burridge-Knopoff model involving stick- 

o ■ 

^v^ ' slip behavior. It is found that cumulative distributions of hypocenter intervals can be described by 

^ ' the g-exponential distributions with g < 1, which is also observed in nature. The statistics depend 

on a friction and stiffness parameters characterizing the model and a threshold of magnitude. The 



Q 

(N 



conjecture which states that qt + Qr ~ 2, where qt and qr are an entropy index of time intervals 
and spatial intervals, respectively, can be reproduced semi-quantitatively. It is concluded that we 
provide a new perspective on the Burridge-Knopoff model which addresses that the model can be 
Q\ • recognized as a realistic one in view of the reproduction of the spatio-temporal interval statistics of 

'^H . earthquakes on the basis of nonextensive statistical mechanics. 

PACS numbers: 05.65.-|-b, 91.30.Px, 05.45.Tp, 89.75.Da 
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I. INTRODUCTION 



Earthquakes occur as a result of a fracture process and a frictional slip of a fault and are categorized into nonlinear 

^ ' complex phenomena in a non-equilibrium open system. In physics, the theory of complex systems and non-equilibrium 

oo ■ 

fvi I systems has not been established, while in seismology the nature of a friction force acting on fault surfaces has not been 

ZX^ , elucidated. Thus, many fundamental problems remain unsolved. On the contrary, statistical features of earthquakes 

I — I , are well recognized as empirical laws, for example the Gutenberg- Richter (GR) law [l|], and the Omori law [2]. As for 

^^ , spatio-temporal intervals between successive earthquakes, unified scaling laws were proposed [3|, Ij]. 

^^ • Inspired by the self-organized criticality (SOC) [5|, many earthquake models have been proposed and then compared 

P*" ■ the statistical features of events with those of earthquakes in nature [6|] . For instance the spring-block model proposed 

/\ ' by Burridge and Knopoff [TJ] (hereafter referred to as the BK model) is well known as a simplified model of fault 

Cu ' systems. Since Carlson and Langer presented a magnitude distribution which was similar to the GR law in one- 



dimensional (ID) BK model [8|, many simulation studies based on this model have been carried out. Recent works of 
this model are focused on the spatio-temporal correlations [9|jll0|i interoccurrence time statistics, statistical properties 
of time intervals between successive earthquakes [ll[ , and a long range interaction I4I , and so forth. Although the 
BK model simplifies complex fault dynamics, the model is useful for discussion of statistical properties of earthquakes 
because the model reproduces the major statistical features of earthquakes, the GR law and interoccurrence time 
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statistics 

Nonextcnsive statistical mechanics proposed by Tsallis [13| and its applications have been paid much attention 
because this statistical mechanics is expected to provide a unified framework for understanding the statistical proper- 
ties of spatio-temporal correlated systems and complex systems. Nevertheless, earthquakes are phenomena exhibiting 
strong spatio-temporal correlations and great complexity. Hence the statistical properties of earthquakes based on this 
idea have been proposed, such as the modified GR law [ij, lia, [ig , interoccurrence time statistics [17|, ll8|, epicenter 
interval statistics 19|, and scale- free network topology [20|. Especially, Abe and Suzuki found that the cumulative 
distribution of interoccurrence times and epicenter intervals can be described by the g-exponential distribution with 
q > 1 



17[ and with q < I [19|, respectively by analyzing Japan and Southern California earthquake data. Then 



they proposed a conjecture, which states that the sum of an entropy index obtained from the interoccurrence time 
distribution qt and that from th e ep icenter interval distribution q,. is similar to two [19| . This relation is also observed 
by examining Iran earthquakes [l8| . The present author reported that the cumulative distributions of interoccurrence 
times derived from the 2D BK model can be described by the g-exponential distribution with q > I, which reproduces 
the natural seismicity n|. However, as far as we know, spatial interval statistics and the conjecture have not ex- 
tracted yet. It is remarked that distributions of epicenters and network of epicenters in the Olami-Feder-Christensen 



model 



211 were reported 22l . |23| . They showed that these statistical features are similar to the observed behavior of 



real earthquakes [20| . 

We try to work out how the hypocentcr interval statistics arc affected by the change of major physical quantities, 
for instance the friction force acting on the surface of faults, stiffness of a fault, and a threshold of magnitude. To 
achieve our aim statistical properties of the hypocenter intervals between successive earthquakes are investigated by 
analyzing synthetic data created by the 2D BK model. Additionally, we examine the dependence of the statistics on a 
friction and stiffness parameters, and the threshold of magnitude. Finally, the relation between qt and q^ is discussed. 
We can show that the spatial interval statistics derived from the 2D BK model can reproduce semi-qualitativcly that 
statistics in nature. 

II. MODEL 



We display a diagram of the 2D BK model in Fig. [T] In this model, a fault is represented by the 2D network of the 
blocks interconnected by springs, fc|f , fc-f, and kp. In this work, we assume that the slip direction of blocks is restricted 
to the j/-direction. The equation of motion of the block at site {i,j) is expressed by 






fee (y^+ij + y^-ij - 2yjj) + ky{vt,]-i + 2/»j+i - '^■yij) - fc^y*. 



F V 



dt 



(1) 



where m is the mass of the block, y is displacement, and i^ is a dynamical friction force as a function of a slipping 
velocity of a block. In order to rewrite Eq. ([1]) as a non-dimensional form, we define a dimcnsionlcss displacement U, 
a dimensionless dynamical friction force (f), and dinicnsionless time t', respectively, as 

Ui,j = yi,j/Do = ytj/{Fo/kp),F{y,j) == Focfiiyij/vi), t' = ujpt = yjkp/m t. 
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FIG. 1: Block and spring system for the 2D BK model. The model is composed of blocks of equal mass, m, of two different 
of coil springs, k^ and fc^, and of the leaf spring, kp. Each mass is subjected to the friction force, which depends only on the 
velocity of the block. 



where Fq is the maximum friction force and vi is the characteristic velocity at FqJI. Using the parameters defined 
above and Eq. dT]), we obtain 



rf^ 
di'2 



'x(t^«+ij + V,-x.j - 2[/,j) + lt,{U^,J-x + C/,.j+i - 2C/,j) - U^,J - 2a U + 



dV,. 



dt' 



(2) 



where 




V V 

— , 2a=—. 

V* Vi 



■Vp y rvp tOpDo 

Ix and ly are the ratio of stiffness of the a; and y directions, respectively, z/ is a dimensionless loading velocity, a is 
the ratio of the maximum slipping velocity, w*, to the characteristic velocity, vi. We use the non-dimensional friction 
function 4>{U), namely 



0({7) 



-oo,l] 
(1 



C/ = 0, 
[/>0, 



(3) 



{l + 2a[U/{l-a)]} 

where a and a are control parameters, a represents the velocity-weakening tendency and indicates how quickly the 
dynamical friction force decreases with increasing velocity. Simultaneously, 2a is the differential coefficient at [/ = 0. 
In the case for a = 0, the dynamical friction force is constant, 0(0) = 1 — ct, and for large a, the force decreases 
rapidly to 0. According to the relation between the friction force and the slip velocity based on the rock fracture 
experiment [2J], a is on order of unity, a is the difference between the maximum friction force (=1) and the initial 
stage of the dynamical friction force (= '/'(O)). In order to prevent a back slip, which means that blocks slip in the 
—y direction, (j) ranges from — cx) to 1 at [/ = 0. 

Previous works reported that the parameter regime, where the 2D BK model can be recognized as a realistic 
earthquake model, is limited and is estimated to be l-i ~ 1,Z^ = 3, and a « 3.5 llll . 1251 . In this case, the model 

nn rT 

can extract the GR law with 6 = 1 [llj, |25(, a statistical property of stress drop [25|, the constant stress drop, the 



Zipf-Mandelbrot type power law for interoccurrence time statistics 



11 



, and the ratio of the seismic wave velocity 



11 



In this work, the system size is taken to be N^ — 100 and Ny ~ 25. As wc described before, the 2D BK model is 
charaeterized by five parameters, /^,Z^,a, cr, and v. Here the two model parameters, a and v, are fixed to be 0.01, 
whereas 1'^,1't, and a are varied around the optimal parameter regime. It is eonfirmed that hypocenter statistics 
do not change quantitatively by varying a and v. We use the fourth-order Runge-Kutta algorithm with time step 
At = 0.001 to solve Eqs. ([2]) and ([3]) under a free boundary condition. The initial block displacements have small 
irregularities. Earthquake-like events are used after a certain period of time when the initial random configurations 
do not influence statistical properties. A seismic magnitude m in this model is defined as ?n = logj^o ( S"? ^"^i-j ) /1-5, 
where 6uij stands for the total displacement at site (?, j) during an event and n is the number of slipping blocks. 

III. RESULTS AND DISCUSSION 

We study hypocenter interval statistics between successive earthquakes. In this model, a position where a block 
slips for the first time during an event is considered as a hypocenter. The nth hypocenter distance is defined as 
Tn = \T^n+i ~ T^n\, whcrc fn and fn+1 Siie the position vector of the nth and n + 1th earthquake's hypocenter, 
respectively. We introduce the scaled distance, r' ~ r/f, where f is a normalized parameter and is set at 20.0, 
arbitrary. Actually, the statistics we will present later do not change statistically by varying f. 

We focus our attention on the applicability of the o-exponential distribution for the cumulative distribution of 



hypocenter intervals P{> r'). According to the paper 19|, P(> r') can be described by the q-exponential distribution, 
eq{x), namely 

P(> /) = e,{-r'/ro) = [1 + (1 - q){-r' /ro)r^^'-'''>] + , (4) 

where q and rg are respectively called the entropy index and the length-scale parameter, and [a]+ = max [0, a]. It is 
known that the g-exponential distribution converges to the exponential distribution as q —> 1, and that for q > I the 
g-exponential distribution corresponds to the power law [l3|. The inverse function of eq{x) is called the g-logarithmic 
function, hiq{x), given by 

ln,(x) = -^(xi-« - 1). (5) 

1-q 

Then, the hypocenter interval statistics are studied produced by the 2D BK model for different I'^^l'^^a, and the 
threshold of magnitude lUc. In this study, we select the g-exponential distribution with q < 1 for the ideal curve of 
P{> r'), because the probability density function of the hypocenter intervals in this model does not obey the power 
law. 

A. Friction parameter a dependence 

First, we examine the dependence of the hypocenter interval statistics on the friction parameter a. For that purpose 
the stiffness parameters, l'^ and l"^ are fixed [l^. = 1 and l'^ = 3) and a is changed from 1.0 to 5.0. As shown in fig. [21(a) 
and (b), the cumulative distributions of hypocenter intervals are well fitted by the g-exponential distribution with 
q < 1. For a = 2.5 and 3.5, the fitting parameters and the r^-value in other words the correlation coefheient, denoted 
by, i?2 are estimated to be q^ = 0.479, tq = 2.54 and Rl = 0.9988 for (a) and qr = 0.474, tq = 2.56, and R^ = 0.9986 
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FIG. 2: Cumulative distributions of hypocenter intervals between successive events for different a, whereas 1^ = 1 and ly = 3. 
The parameter a is a = 2.5 (a) and 3.5 (b). The makers and the solid lines correspond to the numerical results of P(> r') and 
the ideal curve of P{> r') defined in Eq. Q. We demonstrate that the inset figures represent the semi-g-log plot of the data. 
The fitting parameters qr (o) and ro (x) as a function of a are shown in (c). 

for (b) by means of the least-squire root test, qr decreases as a increases, while tq increases as a increases (see fig. [2] 
(c)). R^ ranges from 0.9985 [a = 4.5) to 0.9991 (a = 1.0). It is found that the hypocenter interval statistics depend 
on a. 

B. Stiffness parameters l^ and ly dependence 

In our second performance of our simulation, the stiffness parameter dependence of hypocenter intervals is studied. 
In this time, a is fixed at 3.5, whereas l'^ and ly are systematically changed. Our results of fitting parameters and 
R^ are listed in Table. HI As can be seen from this table, for large l^ and ly, qr tends to be large, while ro comes 
to be small, indicating the fact that P(> r') follows the g-exponential distribution and depends on the stiffness of 
the system. Thus, we can conclude that the cumulative distributions of hypocenter intervals can be expressed by the 
g-exponential distribution with g < 1 by varying /^, l"^ and a. 

C. Threshold of magnitude rric dependence 



Thirdly, the hypocenter statistics above a certain magnitude nic are examined. We use the data with Z^ = 1, Z^ = 3 
and a — 3.5 because in this case the model yields the realistic magnitude distributions and interoccurrence time 
statistics ll|. For the parameter settings given above the magnitude m ranges from —0.65 to 1.5. To keep the 



statistical quantities accurate, we determine the upper limit of riic to be 1.0. Then hypocenter interval statistics are 
studied by varying rric from —0.6 to 1.0. In fig. [3] we display a change of fitting parameters qr and rg as a function 
of TTic. At first, qr and rg do not change in the magnitude region, —0.60 < rric ~ 0.1. However, for 0.1 ^ nic < 1.0 
as 771c increases, qr decreases and ro increases. In this case the value of R^ is between 0.9924 (for nic — —0.60) and 
0.9986 (for nic = 1-0) so that the ^-exponential distribution is suitable to describe P{> r') in the magnitude domain, 
-0.6 < 77V < 1.0. 



TABLE I: The results of fitting parameters of P{> r') for different l^. and l^ with fixing a = 3.5. 
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FIG. 3: Dependence of g^ and ro on the value of threshold of magnitude rric. The data is employed in the case of l^ 
and a — 3.5. 
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D. System size dependence 



Now, we focus on the system size dependence of the hypocenter interval statistics. We change the number of 
blocks N from 2500 (25 x 100) to 40000 (100 x 400), with fixing l^ = 1,1^ = 3, and a = 3.5; N is taken to be 
N = 2500 (25 X 100), N = 10000 (50 x 200), and N = 40000 (100 x 400) (see Fig. [3]). As shown in fig. H the 
cumulative distribution of hypocenter intervals P{> r') is well fitted by the ^-exponential distribution with q < 1 
for all the cases. The fitting parameters and R^ are estimated to be qr — 0.474, ro = 2.56, and R^ = 0.9986 for 
{N = 2500), qr = 0.474, ro = 5.17, and R^ = 0.9985 for {N = 10000), and qr = 0.477, ro = 10.35, and R^ = 0.9984 for 
[N = 40000). qr do not depend on N , whereas ro increase linearly with N] ro for N = 10000 is about twice as large 
as that for N = 2500, and ro for N = 40000 is about four times as large as that for N ~ 2500. Thus it is concluded 
that the hypocenter interval statistics hold up when the system size increases in the range 2500 ^ A^ ^ 40 000. 
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FIG. 4: The cumulative distribution of hypocenter intervals between successive events for different A'^, whereas l^ = l,ly — 3, 
and a = 3.5. The markers and the solid line correspond to the numerical data and ideal curve of P{> r'). 

TABLE II: Summary of the results of fitting parameters of spatio-temporal intervals between successive events based on the 
2D BK model and on earthquake data. 
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E. Conjecture 



Finally, we check the conjecture demonstrated in Ref. 



3, 



which states that qt + Qr ~ 2. The results of spatio- 



temporal interval statistics by analyzing the synthetic data for different l'^, ly, and a and earthquake data of Southern 

are listed in Table HIl R? in this table means the r^-value of 



California [171 . 



3Q 



Japan 



3Q, 



and Iran 



3 



interoccurrence time statistics. Note that we only focus on the case of control parameters where the (j-exponential 
distribution is more suitable distribution for ideal curve of P{> t) than the exponential distribution (g = 1). In 
addition, for l'^ = l,ly ~ 3, and a = 3.5, the value of qt and R^ are quoted by our previous paper 11|. As shown 
in Table [Til it can be deduced that the sum of qt and qr in this model shows qt + qr ^ 1-5, which is similar to the 
observed value. 



IV. SUMMARY 

We have found a new insight into the 2D BK model in connection with the hypocenter interval statistics between 
successive earthquakes on the basis of the nonextensive statistical mechanics. It is shown that the cumulative distri- 
bution of hypocenter intervals is in agreement with the g-exponential distribution with q < 1, which reproduces the 
observed behavior of real earthquakes. The statistics depend on the dynamical parameters, stiffness of the system, l'^ 
and ly and frictional features of a fault, a, and the threshold of magnitude rric. In detail, the distribution function 
do not change, while the fitting parameters qr and tq depend on /^., ly, a, and rric. Additionally, we have checked the 
conjecture, which is a relation between qt and g,- and revealed that qt + qr ^ 1-5. Although the physical interpretation 
of this relation remains open, our findings can lead to the conclusion that the model is able to extract the spatial 
interval statistics in nature semi-qualitatively. 

It is known that the BK model is highly simplified, so that many features of the fault system have been neglected. 
However, the model can be recognized as a realistic model in view of the reproduction of statistical features of 
spatio-temporal intervals between successive earthquakes based on nonextensive statistical mechanics. In a future 
publication, we will focus on the network of epicenters in this model. Finally, we hope that this work is a first step 
toward understanding fully the origin of statistical properties of earthquakes and of other physical systems exhibiting 
the stick-slip motion. 
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